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Heavy-Traffic Analysis of Multi-Type Queueing Under Probabilistically Load- 

Preferential Service Order 

D. P. Gaver 

Naval Postgraduate School 
Monterey, California 93943 

J. A. Morrison 

AT&T Bell Laboratories 
Murray Hill. New Jersey 07974 

ABSTP\CT 

A model of queueing for a single server by several types of customers (messages, or 
jobs), with a simple dynamic priority rule, is considered. The rule is equivalent to 
selecting the next server occupant type with a probability proportional to the number of 
that type enqueued. The situation studied here occurs in fields such as computer and 
communication system performance analysis, in operational analysis of logistics systems, 
and in the repair of dements of a manufacturing system. It is assumed that the population 
sizes of the items of different types are large, and that the mean service rates are 
correspondingly large, in comparison with the service demand rates. Moreover, it is 
assumed that the system is in heavy traffic. Under these assumptions, asymptotic 
approximations are derived for the steady-state means and covariances of the number of 
items of different types either waiting or being served. Numerical comparisons with 
simulated results show excellent agreement. 
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Heavy-Traffic Analysis of Multi-Type Queueing Under Probabilistically Load- 

Preferential Service Order 

D. P. Carer 
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1. INTRODUCTION 

This paper is concerned with the analysis of queueing for a single server by several 
types of customers (messages, or jobs) under heavy traffic. A simple dynamic priority 
rule is presumed to schedule next service among waiting customers when the server 
becomes available. The rule is equivalent to selecting the next server occupant type with a 
probability proportional to the number of that type enqueued. For generalizations, with 
special attention to logistics and repair problems, see Gaver, Jacobs and Pilnick [1]. 
Problems similar to that considered here have been studied by Towsley [2], Yao and 
Buzacott [3] and by others. 

The situation studied here occurs in fields such as computer and communication system 
performance analysis, and also in operational analysis of logistics systems. It is also 
possible to find it in manufacturing, where each of several types of machines suffers 
occasional breakdown and is eligible for repair. The total productivity of the system 
depends on having a sufficient number of each machine type in operation: the present 
scheduling discipline assists in this objective, although others (to be discussed in a 
subsequent paper) are more effective. Note that the current discipline closely resembles a 
processor-sharing scheduling rule which is frequently used for controlling multiprocessors 
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in computer systems. 

In a specific model formulation, there are r types of demand-producing items, and K, is 
the population size of items of type i. If A r f - ( / ) denotes the number of items of type i 
waiting for service, or being served, at time r, then the probability that a new type i 
demand for service is initiated in (t y t+dt) is Nj(t)]dt + o(dt), so that A, is the 

service demand rate for items of type /. The service time of items of type i is 
exponentially distributed with rate v,-, so that the mean service time is 1/v,-; service times 
are all independent. Thus all queue lengths are finite (A',(f) < K,) and a long-run or 
steady-state distribution will always exist, as is true of the simple repairman problem [4], 

A particular example of the probabilistically load-preferential scheduling rules, e.g. 
considered by Gaver, Jacobs and Pilnick [1], is the following. Let N(t + ) = — ), 

A r 2 ( t -P ) A' r (T + )) denote the state of the system at time t + immediately after the 

service of an item is completed. Then, for N(t + ) =*0, an item of type i is selected for 
service with probability 

c,,Y,(t-) 

9/(N(t + )) = — , (1.1) 

2 Cj\>j( 7 + ) 

; = 1 

where cj> 0 for j €(1,2 r). Such a service schedule gives preference, with 

appropriate \\ eights , to those types of items of which there are more waiting for service. 
If one of the queues gets long compared to the others, then it is likely that an item from 
that queue will be selected for service next. Hence this scheduling rule may be regarded 
as a variant of serving the longest queue first. 
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For N(r) # 0, let I(r) denote the type of item which is being served at time r. It is 
clear from the assumptions made that (N(r), I(r), r s 0) is a finite-state-space Markov 
process in continuous time with integer-valued vector state space. In principle, the joint 
probability distributions 

P{S(t) = n, I(t) = i | N(0) = n(0)} = P,(n, r;n(0)), i € ( 1 . 2 r) , (1.2) 

where j £ (1, 2, r), can be found, given the initial condition n(0), by 

r 

solving a system of r FI (Kj + 1) linear differential equations, the Kolmogorov forward 

equations [4]. All probabilistic quantities of interest can be found from such equations, or 
similar backward equations. In practice such solutions involve extensive computing, so it 
is of interest to proceed otherwise. 

To simplify the analysis, we make two basic assumptions. The first is that the system 
is large , i.e. that the population sizes Kj = a oq, where a >> 1 and a,- = 0(1). Since a fast 
server is needed to accommodate a large system, it is also assumed that v,* = a |x,-, where 
jjl / = 0(1). The second assumption is that the system is in hea\y traffic circumstances, 

T 

meaning that I \,a,/|x,> 1. This condition ensures that the server is extremely unlikely 
i = 1 

to be idle. Under these circumstances, Gaver. Jacobs and Pilnick [1] derived a diffusion 
approximation, and obtained numerical results for the time-dependent problem. 

In this paper we consider the steady-state problem under the same circumstances, and 
derive asymptotic approximations for a >> 1 to the means and covariances of the number 
of items of different types in the system, i.e. either waiting or being served. The lowest 
order approximation to the means agrees with that obtained from the diffusion 
approximation [1]. In this paper we also derive the first order correction term to the 
means. Our approximation to the covariances differs from that obtained by the diffusion 
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approximation, when the service rates v, are unequal. However, numerical results indicate 
that the difference is not very large. 

In §2 we formulate the problem, and introduce generating functions. We then 
introduce the scalings corresponding to a large system in heavy traffic, and look for an 
asymptotic expansion in inverse powers of a. The leading term in the expansion is 
obtained in §2, and the first order correction term is derived in §3. Asymptotic 
approximations to the means and covariances of the number of items of different types in 
the system are obtained in §4. Numerical comparisons, which show that the asymptotic 
and simulated results are in excellent agreement, are presented in §5. In the appendix we 
give an alternate derivation of the lowest-order asymptotic approximation to the means 
and covariances, and indicate how we obtain the first order correction to the means by this 
method. We also derive the lowest-order asymptotic approximation to the joint probability 
density function. 

2. AN ASYMPTOTIC ANALYSIS 

Without special boundary conditions the setup described in the previous section is an 
irreducible finite Markov chain and hence possesses a steady-state or long-run solution 

lim P ,(n, r; n( 0 ) ) = p,(n) , i 6 (1, 2, ... , r) , (2.1) 

l - * 

where n = (/i],n 2 n r ). We note that p,(n) = 0 if n,- = 0, and /?,(n) = 0 unless 

OsnsK, n^O, where K = {K j , AN K r ) . The steady-state probability that the system 

is empty is denoted by p(0). For 0^ nsK, n # 0 and n r K, <7,( n) is the probability 
that, when an item departs and leaves the system in state n, then an item of type i goes 
into service. We assume, for the time being, only that <?,(n) = 0 if /:,• = 0. < 7 ,( 11 ) = 1 if 

r 

n } = 0 for j ~ i and 1 < n, < K h and I < 7 ,(n) = 1 for 0 < n < K, n == 0 and n ^ K. 
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We denote by e ( - the vector with components = 5 J; . Then, in a standard manner, we 
obtain 



2 \jKjp(0) = 2 v yPy(«y) , 
y = • >=i 



2 X/fy — ByJ + V, 

L/-» 



P/(e,) = X,W,p(0) + ^ v ; p/e, + e ; ) , 
y = i 



( 2 . 2 ) 



(2.3) 



and, for n, =£ 0, n =?= e, and 0 ^ n < K, 



2 \j(Kj-nj) + Vi 

L/ = i 



P.-(n) 



= 2 l)p,(n-e ; ) + ^(n) £ v^p/n + e,) . 

y' = i y' = i 



(2.4) 



We introduce the generating functions 



«/( x ) = 2 P / C ")*" 1 

0<n< K 



/i r 

•V . 



(2.5) 



and note that u ; (0) = 0. We now multiply equations (2.2), (2.3) and (2.4) by 1, Xj and 
x" 1 • • • x"\ respectively, and sum on i and n. It is found that 









f 


r 


J_ dUj 


b 


p(0) + 2 «/( x ) 

L i ■ 1 J 


- xj 2 - — 

/ = i ^y , 



r V; 



= I — (l-Xj) Uj(\) . (2.6) 



j - 1 



This equation, which holds for general g,(n) with 2 p,(n) = 1, will be useful later. The 

i = 1 

equation is vacuous if x= 1, so that one of the original equations is redundant. The 
normalization condition is 



p(0) + ^ w /0) = 1 • 
i-i 



(2.7) 



We now consider the particular case 
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<7,(n) = “ 



c,rii 



(2.8) 



Z c,n, 
/ = l 



where c, > 0 for / €(1.2 r), and multiply equations (2.3) and (2.4) by and 



v c/n/x'i' ■ ■ ■ x '/ , respectively, and sum on n. After some manipulations it is found that 

/-i 



O 

2 c i x i ~ 

i = i ox < 



j = 1 



= p (0)Ci\iKjXj + c,x, { 2 —1 “ S c,x, 4-iViUj) 

OX. v X; / OX/ 



OX/ Xj J / = 1 OX/ 



r o r r 3w. 

+ 2 c /' r / T“" 2 V*/ * “*/) - v 



/ = 1 



0 X\ ly = 



OX: 



(2.9) 



We are interested in the means and covariances of the number of items of different types 
in the system, namely £(/i ; ) and E(njn k ) - E{iij)E{n k ). We note, from (2.5), that 



Eitij) = 2 T — (!)> £ ( n j n k) = 2 

i = i dx J / = l 



d 2 u, 

dx.dx k 



du; 

(D + b jk — (i) 

OX/ 



( 2 . 10 ) 



We now introduce the scalings corresponding to a large system in heavy traffic, and let 
Kj = ao-j, vj = a[i.j, j t ( 1,2 r); a>>l, (2.11) 

where 

2 > 1 • ( 2 . 12 ) 

j = 1 

Since, from (2.10), we are interested in the behavior of Uj(x) in the neighborhood of 
x = 1 , we also let 

.xj = 1 - £,/a, uj(x) = <J»/|), j €(1,2 r ) . (2.13) 

Then, from (2.6), (2.7) and (2.9), we have 
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p( o) + 2 ^<(0) = 1 • 

i = 1 



(2.14) 



^(4 (0 > + 2 we] + (> - if) 2 |f) - £ Tirp7*/t> ■ «•*» 



and 



M-, 



+ (, ^ d-^/« 

2 C/1 — 7 Cj X My 71 — T~7 

/ = ! ^ a > oij , = 1 (1 



plO) 

a 

1 2 



c,\ 



aj>, ^ (1-^/fl) cji.-ib, 

j'- 

it\ 1 

i- 



a) dtj a(l-£, a) 

(i - -) - 1 £ v>,e, £ c < (i ■ -) lr 

v a ^ a f=i , = j ^ a ) aii 



£, 



1 2 

a ; = : 



k) 






a J 


Hi 




i 


2k) 


d\bj 


i 


a J 





+ - 2 M; (l " ~) 2 c i (l “ “ ) T7TT 
a fix K a ’ / = i V a ' a£;3& 

We assume an asymptotic expansion of the form 

*!>,(€) ~ d4 0) (£) + - + • • • . 



= 0. 



(2.16) 



p(0) ~ 0 , 



(2.17) 



since we expect p(0), the probability that the system is empty, to be exponentially small in 
a. Then, to lowest order, from (2.16) we obtain 

dd4 0) 



' diM 0) ' 

M-i 2 Cl — Ci X M-y 

/ = 1 °?/ j = 1 



= 0. 



(2 18) 



These equations are satisfied by 



4 0) = 



ci <56 



M-. °ii ' 

where the function 6(£) is to be determined. But. from (2.15) and (2.17), we have 

<5d>i 0) 



(2.19) 



£M,£M»--iB-£Mi*r 



;= i <=i 



°kj ' ; = ! 



(2.20) 
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We look for a solution of the form 



0( £) = 9(0) exp (— 2 P/£/) - 

i = i 



( 2 . 21 ) 



so that, from (2. 19), 



vl»i 0 ) (g) = - 



C/B, od;i 0) dB, 

— 0(C). -7T-=— P;9(C) 

M-i °€/ M-i 



( 2 . 22 ) 



We let 



a = 2 



dPi 



/ = 1 M-l 



(2.23) 



Then (2.20) is satisfied if 



A \j(ocj - P;) = C/P;> i-e., P, = 



A \jOL i 

(A \j+Cj) 



(2.24) 



Since we want (3y > 0, j € ( 1 , 2, ... , r) , we want A > 0. But, (2.23) and (2.24) imply that 



F(A) - 2 






; = ] \i.j(A\jTCj) 



1 = 0. 



(2.25) 



But F(0)>0, because of assumption (2.12), and F( + *) = — 1. Since F (A) is a 
decreasing function of A for A > 0, it follows that there is a unique solution A > 0 of 
F (A) = 0. From the normalization condition (2.14), and (2.17), we have 



2 ^ 0) (0) = i 

i=i 



(2.26) 



It follows from (2.22) and (2.23) that 



0 ( 0 ) = - 



(2.27) 



This completes the determination of 0((j), and hence vhf 0) (^). 
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In the next section we determine il/^(£). The reader who is not interested in the 
details may proceed to §4, where the asymptotic approximations to the means and 
covariances are evaluated. 

3* THE CORRECTION TERM 

We now consider the first order correction term in the asymptotic expansion 

(2.17). It follows from (2.16) that 



4-/ 2 c i 77 “ c i 2 ^y 



/ = l 



y=i 



ovl/^ 1 



r 

M-i 2 c &i " : , + c i S M-yC^y - ^) 



/ = i 






C/fi/il/} 



( 0 ) 



2 2 ci 

}- i '=i 

od/} 0) 



;' = i 
od; (0) 

-T- - 2 oWi 01 

y = i 

(3 2 vl;; 0) 



2 C A/ 2 Mf 2 c l -c 



; = i 



- J J = i / 



= i 



It is found, from (2.21), (2.22) and (2.24), after some reduction, tha' 



O) 



^ oUJ '" ^ 

^ 2 ci — ci 2 Mv -77- 

/= 1 i = l 



C/Pi 

p./ 



0(0 



; = i 



X CJ P;-Ac, (J.,- + A\X, 2) c ; P;(2£, — £,) - 2 2 f/PyC; 

V=1 y= i /=1 ; = >• 



. (3.2) 



We note that if we sum (3.2) on i , and use (2.23), then both sides of the equation are 
identically zero. 

In view of (2.21), we let 



0^1 n (0 = xl n (0 ex pf“ 2 Ptdt) 

A = 1 



(3.3) 



and we define 



fi=2 <v (3 .. D = 2 c) fJ. 



( ? 4 ) 
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Then, since 9(0) = — 1M, 



c,P. 



,(1) 



A-\±j 



D A c ,• p. , A5p./£, ~ (2A)J., S) ^ 

7-1 



(3.5) 



We look for a solution of the form 



xJ 2) (€) = - 

P-. 



o6 (1) 

• of,- 



- P,6 (!) (€) 



a, + 2 b * 



(3.6) 



*«i 



and note that the terms involving d> (1) are annihilated by the operator on the left side of 
(3.5), i.e. they provide a solution of the homogeneous equation. 



It is found that (3.5) is satisfied by (3.6) if 



C/P/ 



cfii 2 V-jbjk ~ B\Xjb ik = —7— [(fi -2Ap./)c*P* + ABp.,5,*], 
j = l A*p, 



(3.7) 



and 



c iP. 2 v-j a j - B = c,- 2 M-y b ji - M-/ 2 c i b ii ~ ~f' ' (Ac ( -p.,--D) . (3.8) 



;=i 



; = i 



/ = l 



AV; 



If we sum (3.7) and (3.8) on /, and use (2.23) and (3.4), then both sides of each equation 
are identically zero. Since we are looking for a particular solution of the inhomogeneous 
equation, we may take 



2 V-jhjk = 2 V-pj = 0 



(3.9) 



; = i 



;' = ! 



Then (3.7) and (3.8) imply that 
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and 



bik ~ 



a, = 



c.P, 

AB p. 



-[(2—2-1 

, n Ap.,) 



Ck$k ~ 58, 



2c, p, 
AB 2 ^ 



( D —Bcj) . 



It remains to determine the single function 6 (1) (£) in (3.6). But, from (2.15), 
and (2.20), we have 



r r ( IV dUT; \ r 

2 M; 2 + ~T~) ~ 2 



V V 
;=i /=1 



^ ' , = 1 






- 2 M; 2 -rr- * 2 MW • 



j-i i-i ;-i 



It follows, from (2.21)-(2.24), (2.27) and (3.3), that 



,<i> 



1 r r r r dV s 

x 2 2 x!» + 2 Mj 2 - 2 »*1" 

^ ; = 1 i = l j = 1 i = l 



^ ; = i 



T 2 

A ; = 1 



Then, from (3.6). we obtain 



Vw vi e 2 <t> (1) , i ^ c - a<f> (1) ' . . ,, 66 (1) 

2 M; 2 ae ae + a 2 C ./PA/ 2 2 (AX y +c ;^; ^ 

) = 1 1 = 1 l 1 ' A j = 1 1 = 1 r-1 °S| _/ = 1 °?7 



= 7 2 ( c ,- a \ ; )^ ; 2 -2 2 C ; 3 ; e ; i ( o ,+ i ^ 1*1 

^ j = 1 A ; = 1 1 = 1 *-l 7 

“ 2 M; 2 bij + 2 M- ; €/( a ; ~ 2 bjktk) ■ 



j = 1 1=1 j = 1 



* = ] 



We look for a solution of the form 



6 (1) (|) = k -r 2 'Ytit T { S i Ukitkti 

k = 1 “ * = 1 /=1 



(3.10) 

(3.11) 
(2.17) 

(3.12) 

(3.13) 

( 3 . 14 ) 

(3.15) 



where we assume, without loss of generality, that 
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^kl ~ <«> Ik ■ 



We define 



and 



gjk 



Cj$, 



2 bik V-jbjk + 4 



i*l 



r 

T; = S — “y. 

, = l M-/ 



hj - + 7 , ; 2 a i + 'S' 2 bij V-jOj , 

^ / = i i' = i 



«- 2 

i = i 



c,7/ 



M-i 



Then, it follows from (3. 14)-( 3.20) that 



2 2 

j- 1 A=1 



C/Py 



r* - (AKj+Cj)oi jk + 



2 ij 

7 = 1 



^ c yP/ 04 



’ c j)y. j + h j 



In view of the symmetry in (3.16), we deduce from (3.21) that 

[(AXj + Cy) + (AX k + C k )]d>j k = -^(cfijTk + CtPify) + g jk 4- g kj . 



Hence. 



(c ; (3 y r t + c k ^ k Tj) + 

Jk A[(AX^ + c ; ) + (AXjj+q)] Jk 

where 



, _ (gjk + §ki ) _ 

Jjk [(AX, + c y ) + (AX*i-c*)] Jkj 



(3.16) 

(3.17) 

(3.18) 

(3.19) 

(3.20) 

(3.21) 
= 0. 

(3.22) 

(3.23) 

(3.24) 



We define 



- 13 - 



c = 2 

1=1 M-r 

Then, from (2.23), (3.4), (3.10), (3.17) and (3.24), we obtain 

gjk = yr cfijc&k (- L - y) + , 

A- J K 4, M-Jt 

and 

3( (^y(3 jbjK — ^yPy^tP/t) 

[(A\j + Cj) + (A\ k + c k )] 



(3.25) 



(3.26) 



(3.27) 



From (3.18) and (3.23) we have 

r = JL y c (c,(3.r t + cip^r,) + y fj_ , n~>$) 

k A ft \Lj[(A\j + Cj) + (A\ k + c k )] ft, V-j Jk ' 

which is a linear system of equations for T*., k € (1 , 2, ... , r), with known coefficients, 
which has to be solved numerically, in general. Once the T*. have been calculated, c oj k is 
given by (3.23). Also, hj is given by (3.19). We define 

l 

G = T (3.29) 

,- = l F, 



Then, from (2.23), (3.4), (3.10) and (3.11), we obtain 



hj XjTj 




PjMCy-Oj. 



Next, from (3.21), we deduce that 



(3.30) 



(AXj + Cjjyj = j CjPj + hj. 



(3.31) 



Hence, from (3.20), we find that 



. U <7P.» 

A 



A , UjiAXj-^Cj) 



Cjhj 



V-jtAXj + Cj) 



(3.32) 
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But, from (2.23), 



1 - T 2 



fa 



= X 



c Jhhi 



> o. 



(3.33) 



A “ i \ Lj(AKj + Cj) j = i \ ij(AKj + Cj) 

Hence e is determined, and then 7 j is given by (3.31). 

It remains to determine the constant k in (3.15). But from the normalization condition 
(2.14), and (2.17) and (3.3), 



0=2 vi»i !) (0) = 2 xl n (0). 

1=1 1=1 



(3.34) 



It follows from (3.6) and (3.15) that 



2 

i = ! 



— (y,— P/K) + a, 



M-i 



= 0. 



(3.35) 



Hence, from (2.23), (3.11), (3.20) and (3.29), we obtain 



k = — ( AD - BG ). 

^ D - 



A A 2 B 



(3.36) 



This completes the determination of 6 (1 '(£). and hence X/ 1} (€) 2 nd th| ]) (^). 

4. THE MEANS AND COVARIANCES 

We now evaluate the asymptotic approximations to the means and covariances of the 
number of items of different types in the system. From (2.10) and (2.13) we have 



dll), ^ r _ / x ^ d 2 ll); 



Eirij) = - a 2 T~ (0). E(njn k ) = a 2 2 TZ ' (0) - &/*£(«;). (4.1) 



1 = 1 °tj 






Hence, from (2.17), we obtain the asymptotic approximations 



EOlj) - - 2 
( = 1 L 



oil); 0) 5il;l n 

a -tt-( 0) + -rf-(O) 









(4.2) 



and 
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E(njn k ) - b jk E(nj) - ^ 

1 = 1 



, <5 2 ii>| 0) a 2 ^ 15 

a- ~ tz (0) + a — - 1 (0) + 



SBjdh 



HjHk 



(4.3) 



From (2.22), (2.23) and (2.27) it follows that 



r ao»i 0) 

- 2 — (0) = (J, 

i = 1 “V 



(4.4) 



Also, from (3.3), (3.6) and (3.15), we have 



oili ^ ^ C ' 

' (0) = — (C 0 y-p, 7> ) - bjj - p ; x, (1) (0). 



<35, 



M-I 



(4.5) 



Hence, from (2.23), (3.10), (3.18), (3.25) and (3.34), we obtain 






i=i °$j 



(0) - r, - Ayj + 



(4.6) 



From (4.2), (4.4) and (4.6) we find the asymptotic approximations to the means. 



E(rij) ~ a(3, + 



A yr r i +C fii (7^ + fr - Jj 



l)} 



(4.7) 



The quantities |3y are given by (2.24), where A is the positive root of (2.25), and B and C 
are given by (3.4) and (3.25). Also Ty, j € (1, 2, .... r), satisfy the linear system of 
equations (3.28), subject to (3.27), and 7 y is determined from (3.29)-(3.32). 



Next, from (2.19), (2.21), (2.23) and (2.27), it follows that 



r <3 2 *1/J 0) 

2 iTTr< 0 > - M 3 *' 



, = 1 dtjdtk 



(4,8) 



Hence, from (4.2)-(4.4) and (4.8), we obtain 



E(n.n k ) — E(n;) E(n k ) 






a I P; ~ 2 

l i = l 



2, 1.(1) 



o-u,; 



o£y 3£; ; 



7^(0)- (3 






in 









ail; 

(0) + p*— — (0) 

otj 



(4.9) 
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But, from (3.3) and (3.34), 



2 
i = l 



„ ao/< ]) /rt . 

T7~b~ ( 0 ) + 3, ~T~ <°) + P* ~ — (0) 
dEjotk °s* 



otj 



= 2 



' d 2 x\ h 



= 1 ofy bt k 



( 0 ) 



(4.10) 



Also, from (2.23), (3.6), (3.15) and (3.16), we have 



' a 2 xl n 

2 ^rV' 0 ) = 

i = l °lj°ik 



= — A to 






(4.11) 



From (4.9)-(4.11) we find the asymptotic approximations to the covariances, 



E(njn k ) - E(tij)E(n k ) ~ a($j 5 Jk - Ao> jk ) + • • ■ , (4.12) 

where c o jk is given by (3.23) and (3.27). We note that the covariances, as well as the 
means, are of order a. 



We now consider those systems for which 

Xj = X. Cj = c, ; €(1,2 r) , (4.13) 

since it is then possible to explicitly solve (2.25) for A and (3.28) for T*. We note that 
assumption (2.12) is now 



a,- 



x 2 — > I- 

7 = 1 ^ 



The solution of (2.25) is found to be 



(4.14) 



Also, from (2.24), 




AX a, 

3 , = — , 

J (AX-rc) 



(4.15) 



(4.16) 



and, from (3.27), 
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fjk ~ 



(AX -He) 



K - 7 f &>) 



where, from (3.25), 



„ ^ P« 

C = c ^ — 

i = l M-f 



We define 



r r. 

n = 2 — 

j = \ v-j 



Then, from (3.28) and (4.15)-(4.17), 



A(2AK + c)T k = c : n(3|. + 2cA$ k (— - -^■1 

y P-k A‘ ' 



If we divide (4.20) by p* and sum. we find that 



fl = (AX — c) 

A 2 X 



Hence, from (4.20), 



r t = 






(2AX-rc) 



2X cC 



M-t A J X 



(AX + c) 



It follows from (3.23), (4.17) and (4.22) that 



P / b jk - A(X>j k 

cp/5,|. 






(AX — c) (AX + c)(2AX + c) 



C 



A 3 X 



(2A : X 2 t2AXc-c 2 ) - 



( 1 


1 Xl 


— - — 







which gives an explicit expression for the asymptotic approximations (4.1 
covariances. 



(4.17) 

(4.18) 

(4.19) 

(4.20) 

(4.21) 

(4.22) 

(4.23) 
2) to the 



It remains to calculate the first order correction term to the means in (4.7). The 
quarries h Jt e and ~/j are obtained in a straightforward manner from (3.29)4 3 *3) and 
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(4.22). It is then found from (4.7) that 

c3 Py ( 1 C\ 

E(n ^ ~ + A(AX + c)(2AX-c) I vj aJ" 



(4.24) 



Since our approximation (4.12) to the covariances differs from that obtained by the 
diffusion approximation [1], we give an alternate derivation of our results in the appendix. 
Corresponding to the scalings in (2.11) we let iij = a{ij -r \Ta v y in (2.4), with <?,(n) given 
by (2.8). We then develop an asymptotic expansion in inverse powers of \fa for 
<b,(v) = a rl2 pi( n). It is found that, to lowest order, p,,-6,-(v) ~ c,(3, 0 (0) (v), where 
A <t> ,0) (v) is a multivariate Gaussian probability density function. From this we are able to 
obtain the lowest order approximations to the means and covariances, and we again obtain 
(4.12) and E (nA = a(3 y + 0(1). Although we omit the rather lengthy details, we have in 
fact carried out the analysis to the next order in the asymptotic expansion, and have 
verified that it leads to the approximation (4.7) to £(/i y ). It is of interest that the 
asymptotic expansions of the densities are in inverse powers of "S fa, whereas the 
expansions of the moments are in inverse powers of a. 

The difference in the equation for the covariances between our approximation and the 
diffusion approximation [1] is pointed out in the appendix. 



5. NUMERICAL EXAMPLES 

In order to check the accuracy of the approximations proposed, several systems were 
both simulated and approximated. In all cases studied the agreement between the 
asymptotic expansion approximation and simulation was excellent for both mean and 
standard deviation of the steady-state distribution. The mean of the diffusion 
approximation agrees with the lowest-order asymptotic approximation; the standard 
deviation of the current diffusion approximation agrees precisely with the lowest-order 
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asymptotic approximation only when the service rates are identical, but correspondence is 
quite close numerically for all situations investigated to date. 



Here are the systems, and their properties. 



i: 1 

K t \ 100 

X, : 1.3 

v,- : 50 

c,-: 1 



System 1 
(r — 5) 

2 3 

110 120 

1.3 1.3 

100 300 

1 1 



4 5 

130 140 

1.3 1.3 

400 450 

1 1 



Means 

Asymptotic: 
Simulation: 
(95% Conf.) 

Std. Deviations 

Asymptotic: 
Diffusion: 
Simulation: 
(95% Conf.) 



81.40 

81.14 

(81.04, 81.24) 



3.73 

3.75 

3.75 

(3.67, 3.82) 



89.54 

89.39 

(89.23, 89.54) 



4.38 

4.38 

4.38 

(4.31, 4.45) 



97.68 

97.47 

(97.37, 97.57) 



4.93 

4.91 

4.89 

(4.80, 4.97) 



105.82 

105.60 

(105.42, 105.79) 



5.23 
5.20 

5.24 

(5.17, 5.31) 



113.96 

113.74 

(113.53. 113.95) 



5.50 

5.47 

5.56 

(5.48, 5.65) 



System 2 
(r = 5) 



i : 


1 


2 


3 


4 


5 


K r . 


50 


100 


150 


200 


250 


k,: 


1.3 


1.3 


1.3 


1.3 


1.3 


v,-: 


50 


100 


300 


400 


450 


c, : 


1 


1 


1 


1 


1 



Means 


Asymptotic: 


39.18 


78.37 


117.55 


156.73 


195.91 


Simulation: 


39.10 


78.17 


117.35 


156.26 


195.57 


(95% Conf.): 


(39.03, 39.18) 


(78.05, 78.28) 


(117.12, 117.57) 


(155.95, 156.57) 


(195.30, 195.8 


Std. Deviations 


Asymptotic: 


2.78 


4.25 


5.85 


7.13 


8.34 


Diffusion: 


2.85 


4.33 


5.90 


7.24 


8.56 


Simulation: 


2.78 


4.21 


5.90 


7.19 


8.33 


(95% Conf.): 


(2.75, 2.80) 


(4.16, 4.27) 


(5.81, 6.00) 


(7.13, 7.25) 


(8.26, 8.40) 



The simulations were carried out on an IBM 3033 computer at the Naval Postgraduate 
School, using the LLRANDOMH random number operating package; Lewis and Uribe 
[5]. Time-dependent queue lengths were simulated: an event clock was advanced at either 
job arrivals or service completions, and queue lengths were suitably incremented or 
decremented. The current queue lengths were recorded at fixed discrete time steps; 500 
independent replications were recorded. A batch mean process was utilized to obtain the 
confidence limits. For further details see Pilnick [6], 

For the two systems displayed, and for many others explored, the leading term 
asymptotic agreement with simulation was well within the 95% confidence limits displayed 
for the latter. Diffusion approximation was also good, but slightly less accurate than the 
asymptotics described here. A full description of the diffusion approximation approach 
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taken is given in Gaver, Jacobs and Pilnick [1]. The small size of the standard deviations 
as compared to their means is noticeable: certainly there is little resemblance of the current 
system behavior to that of a system with multi-type Poisson arrivals for service. In the 
latter situation the queue length standard deviation will be close to the corresponding 
mean queue lengths in heavy traffic. Of course in the case of multi-type Poisson arrivals 
there will be a steady-state solution only if a suitable traffic intensity for the system is less 
than unity; no such condition need be satisfied for the finite systems considered in this 
paper. 
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APPENDIX 



We here give an alternate derivation of the asymptotic approximations to the means 
and covariances of the number of items of different types in the system. Corresponding to 
the scalings in (2.11), we let 

n j = a + Va \’j, pj( n) = $> i (\)/a r 2 , ( A 1 ) 

and it will be shown that (3 is given by (2.24). From (2.4) and (2.8), we have 



nA-(v) - i \j(ccj-£j-^L) [<h,(v) - <b'(v-ej/Va)] 

j = l v x a J 

1 A /— Cj($t + Vi/\/~ci) JL r- 

- — 2 ^'6,(v-e y /Va) - — 2 M-j*;(v-re ; /\ a). (A2) 






2 c t (pit + v t /Va) 1 1 
A' = 1 



If we expand in inverse powers of \ a . we obtain 



r ( \’i \ r i d<bj 

„,*,(>•) - 2 («;- 0 - [— — (v > + 



£/_ 

B 






V a 



B k % 



2 P, 

; = i 



- — 2 M^'( v ) + 

G ; = i 



1 06, 

*,(>') + — — -r— (v) + 
V a o\’j 



, (A3) 



where B is given by (3.4). 



If we sum (A2) with respect to i, we find that 

2 P/<5>x(v) + 2 =^1 2 [cb,(v)-<t>,(v-e ; -/Va)] 

/ = l ; = i ’ k ; . = i 

= ~ 2 x j 2 6,(v-e ; 'Vfl) - 2 M-ycb y (v-e,'V2) . (A4) 

° ;=: ' 1=1 ;=1 



It follows that 



A-2 



r , V, \ r r 56, 1 

2 vr) 1“ v) ~ 2 V? 



1 o 2 6i 



rr(v) 



OV- 



2 X; 2 M>i(v) 



Va ,«1 v , = 1 



] + 2 m-7 

7 = 1 



56; 1 5 2 6; 

^(v) + -y — ^(v)- • • • 

ov i 2v a av~ 



We assume an asymptotic expansion of the form 



6,(v) ~ 6^ 0) (v) + —4 — <J>| n (v) + — 6( 2) ( v ) 
V a a 



Then, from (A3), we obtain 



ix,-6( 0) (v) = C ~y- 2 M-;6j 0) (v) = c,(3,<t> (0) (v) , 






where 4> (0) (v) is to be determined. But, from (A5), 



r r r <36/°^ 

2 M«/-W 2 — = 2 -r 

j — \ i = l °'J 7 = 1 0 7 



It follows from (2.23) and (A7) that 



' 6$ (0) 

2 MX,(o,,-fS,)-f,|3,] -r— - 0. 

j = l 0 V 



This equation is satisfied if (2.24) holds, which we assume to be the case. 



Next, from (A3) and (A6), we have 



63, 



M>J n (v) - 2 M-;6j n (v) = ~ fv, — 2 4t v *l 2 M- y 6j 0) (v) 

D j=i a 3 u = i ;=i 

c,(3, ' 56} 0) ' 56$ 01 

+ — 2 ^ ~^r ~ 2 x / a 7-P;) 

D 7 = 1 0 7 ; = i °v 



• (A5) 



(A6) 



(A7) 



(A8) 



(A9) 



(A10) 



Hence, from (2.24) and (A7), we obtain 



(All) 



n.<f>i n (v) - c,(3,<t> <n (v) + ci>>)o<»(v) 






dO (0) 

J * v j ' 



where is to be determined. But, from (A5) and (A6), 



2 2 



' 66^ ]) 

“ 2 M-y 



7 = 1 



, = i dv ; j = i 



3v, 



r _L d&\ 0) j r 

= 2 ^j v j 2 7~ r ~ 2 2 



5 ; 6< 01 



7 = 1 i = l ^ v 7 



7 = 1 



f = 1 OVj 



r r y r 3 ^ d) 

+ 2 Xy 2 ^V) + 72 *J— T-- 



7=1 '=1 



l = 1 ov- 



(A12) 



From (A7) , (All) and (A12), with the help of (2.23), (2.24) and (3.25), it is found, 
after considerable reduction, that 




- t 2 — 2 cfij 4 -( v k < p (0) ) - ± u \ r cj ) 4 r ('- j & 0, ) = o. (a i 3 ) 

A k = ] M-* 7 = 1 0X j 7 = 1 ° X j 

an equation involving <f> <0) only. From (Al) the normalization condition implies 
asymptotically, from the Euler-Maclaurin summation formula [7], that 



X X 

f • • • / 2 *«('■)<*’ 1 • • 

- X - X / = 1 


X 

• dv r = / 2 <±>/( v )<^ v ~ 1 ■ 

-X I = 1 


(A 14) 


Hence, from (A6). 






X 

f 2 6 \°Hx)dv = 1, 

- X 1=1 


X 

f 2 <b| n (v)rfv = 0. 

- X ; = ] 


(A 15) 



It follows from (2.23) and (Al) that 



A-4 



f cp(°)(v)^v = — 

a 



(A 16) 



We let 



x 



Mi = A f v,Q (0) (v)dv. 



(A 17) 



Then from (A 13), after multiplication by v,- and some integrations by parts, we obtain 



It follows from (3.33) that A = 0, and then from (A18) that M,- = 0. Consequently, from 
(Al), (A6), (A7), (A15) and (A17), E(nj ) = a(3y + 0(1). We have in fact carried out the 
analysis to the next order in the expansion (A6), and have verified that it leads to the 
asymptotic approximation (4.7) to E(>ij). However, we omit the rather lengthy details. 

Next we consider the covariances, and let 



(A\i + Ci)Mi = Ci \ 3, 4- A =£ — M k . 



k - 1 V-k 



(A 18) 



Hence, 




(A 19) 



M a = A f \'i\’,& 0) (x)dv. 



(A20) 



If we integrate by parts and use (A 16), we obtain 




(A21) 



and 



x 




(A22) 



Hc^'e. from (A 13), we find that 



A-5 



[(AX, + c,) + (A\i + Ci)]Mjt - -j- ^ — (CjpjM/jt + c/P/A/,^.) 

^ i=i M** 



= 2c,(3,S;/ 4- — Cj 

A 




2C 1 1_ 

A M-i M-; 



)■ 



(A23) 



If we let 



Mu ~ P/ 5,v — A to// , 



(A24) 



then we obtain equation (3.23) for a where T ; and f jk are given by (3.18) and (3.27). 
The asymptotic approximations (4.12) to the covariances follow from (Al), (A6), (A7), 
(A20) and (A24). 

We remark that in the diffusion approximation [1] to the covariances, the off-diagonal 
terms on the right-hand side of equation (A23) are absent, although the diagonal terms 
agree exactly. 

We will now show that AO ,0, (v) is a zero-mean multivariate Gaussian probability 
density function. We introduce the characteristic function 



x <°)(y) = a / e ix 5 O (0) (v) d\, v y = 2 l 7>7 - 



(A25) 




where i = V— 1 . If we integrate by parts, we obtain 




(A26) 



and 




(A27) 



It follows from (A 13) that 



x ( 0 ) (y) 



2 c j Vjyj “ T 2 c ./ 3y(“" “ t) 2 3* W 

L j = i A j-i K ^J * ' *«i 

1 r Qt r r 

- 7 2 7T 2 ^,y, %- + X M^+<v>% - 0. (A28) 

A *-l M-* 7 -i d y* ; = 1 d >7 

But, from (A16), ( A 17) . (A20) and (A25), since Mj— 0, we have 

ov (0) a 2 v (0) 

X (<3, (0) = 1, — (0) = 0, - t-V-(O) = M Jk . (A29) 

dyj oyj ay k 



It is straightforward to verify, with the help of (A23), that (A28) and (A29) are satisfied 
by 



X (0) (y) = exp f- zr X 2 Wfc,y,yJ , 

Z / = 1 m- 1 



(A30) 



which establishes the desired result [8]. 
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